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SUBJECT: PROGRAMS FOR SOLVING SECULAR EQUATIONS 

To J Scientific and Engineering Computation Group 

From* FoJ Corbato 

Introduction 

Production programs for Whirlwind I are now available for 
solving two types of secular equations by use of the techniques described 
in a report by A Jfeckler and are summarized briefly below. These are 
the ordinary variety , (case l) 9 ^_ ^"^ik * Yik^k* and tile & eneral 

variety, (case Il)$>Z_ %•* Yik * ^— S ii^ik^k' where k * 1 » 2 > .•o,'n, 

H and S are real symmetric matrices and S is positive definite- The two 
programs 9 for case I and for case II, can handle matrices of order 
1 4 n 4 32 o The results are given ^photographically and consist of the 
input data and the A, and W., where in case I, 2 Yi1 rik * °ik 



J 



and in case II, 2 ^/i S iiYiv s ° 0v° *n addition the intermediate 

■i» J 

results of case II 9 if desired may be displayed. 



Ifethods of Solution 

Matrix diagonalization is the elementary process used to solve 
both cases. The diagonalization of a symmetric matrix M-, ■ M, - is 
accomplished by successive 2 by 2 "rotations" of all the matrix elements 
associated with the indices i and j where M.. is the largest off -diagonal 
element o It can be simply shown that such a process converges « The ex- 
plicit transformations applied ares 

°A«» Meckler 9 Quarterly Progress Report 9 Solid-State and Molecular Theory 
Group, MoIoTo, Octo 15, 1954, P° 15. 
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M /k - M * k 



M ik ■ cM ik + s V 



M 






" sM ik + oM jk 



k ^ i,j 



M ii * C ~*L1 + 2soM ij + 8 \ 



M 



JJ 



"^ii " 2SCM ij + °\ 



M iJ ■ ° 



where 



^ll + t 2 



te 



and 



t - / < M ii - M ii> * N.«li - M i.i) * *<*ur >W 



»2M 



ij 



+2M 



11 



V 



(M i± - M^) . Vjtt^ - *, )* ♦ 4(M ±j )* f (M^M^) 



Similarly , the unitary transformation, U , (initially a diagonal 
unit matrix) , is modified after each 2 by 2 rotation by the corresponding 
transformation affecting the ith and jth columns. 



U 



km 



U 



km 



m 



/ i.j 



V - cU ki* sn kj 



% = - str ki + oU kj 



This process continues until M.. ^ 2 , where c is the preset criterion 
value, the result being that the diagonal of the final M. . (i.e. the eigen- 
values , X ) and the final U satisfy the conditions : 



w 



y^ M o,Uo, - U 4u^t, or ^ = u ^ (where A does not 
j=l ^ JJC ^ K multiply as a vector 

but as a scalar) 
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In nearly all applications* the original matrix, M, was computed 
from an operator,'^ , by means of a set of basis states « These states 
can be represented as a set of unit vectors, e. ( in a Hilbert space, and 
what is desired is the transformation which when applied to the original 
basis states gives a new set of basis states that when applied to the 
operator nv , yield a diagonal matrix,. To obtain this transformation it 
is noted that if 



and 



M ij - <^M^> 



i ^K c 2_ ijJm. .u,, 

°£k k J-j jfi ij jk 



then % A k W ^iVi 1 ' j/ F ik 






<e£M"0 



Zut 



Thus the desired transformation is e, m /__ U, .e.» Stated in another 

k j~ k ^ ^ 

manner, the kth eigenvalue of H is associated with an eigenvector consisting 

of the kth column-vector of U., where the components of the eigenvector 

refer to the original basis states used to compute M„ 

For case II, the procedure of solution is the following First 

t t -1/2 

S is diagonalized so that SUg,» tJg,S , where S is diagonal. Then S ' 

is formed where S™ ' 2 « U gf (SO" 1 ' 2 !^, and the new matrix H § is formed by 

matrix multiplication where H ■• S"" ' TIS ' • Then a second diagonalization 
is made such that 

and finally Y * S^Ajg, where it is seen that Y f S Y " *• 

To verify that this A and Y form desired solutions, it is 
observed that 
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ss'-^^.X 



STX 



Specification of Input Data 

Both programs have their input data, (also called parameter 
sets) 9 located in the same storage locations and moreover do not destroy 
during operation any of the necessary input information. Thus a series 
of secular equations each differing from the previous by the addition of 
a row and column could be conveniently solved by supplying the complete 
matrix (or pair of matrices in case II) in the first parameter tape and 
then merely changing the matrix order 9 n 9 for each subsequent parameter*, 
Furthermore 9 a parameter tape for case II could be used with program I 9 
but would 9 of course 9 give solutions for case I» (Program I would also 
disturb the storage of the unnecessary S matrixo) 

In the listing of the input data 9 two types of number conventions 
are usedo These are the single=register (15 9 0) integer , (less than 32768 
in magnitude) , which may have a sign but no decimal point 9 and the double- 
register (24s>6) generalized decimal number which must have both a sign 
and a decimal point* The specific locations of the input data ares 



34| 


+n 


35| 


+c s 


36| 


+C H 


37| 


+k 


l4 


+H 11 



matrix order 9 l^n^32 

diagonalization criterion of S 9 Cg^-60 

diagonalization criterion of H, c H >-60 

identification number _> 

listing of the H matrix | 

(continued on next page) 



. (15*0) 
integers 



DCL-58 



3104 



+H 



12 



+H 



22 



+H 



+H 



13 

23 



+H 



nn 



+S 



11 



+S 



12 



+S 



22 



♦S 



+S 



13 
23 



+S 



nn 
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listing of the S matrix 



(24,6) 
numbers 



Both programs start at register 32 decimal and stop on an "si n 
instruction in register 33 decimal (41 octal). If many sets of parameters 
are to be included in one big tape and are to be automatically run con- 
secutively, all but the last set should contain an additional 33| sp26, 
and the last set should contain a 33| si0» Each set on a multiple- 
parameter tape must be preceded by a tt fc H and then a n (24»6) tt § the per- 
formance request should have, under the heading of operation instructions 9 
a '♦Turn off si 1 switch** on the line preceding the multiple-parameter 
tape operation instruction « 



Criterion Values 



The values of C„ and G g are used to terminate the diagonalization 

procedures after sufficient accuracy has been obtained » For both the H 

and the S diagonalization 9 each c should be chosen such that 

Imaximum absolute error I i maximum absolute error! 



K-\ 



j Imax 



• i in eigenvector components! ^ I in eigenvalues 
* X l>klaa* 
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The accuracy of the computations is such that the right-hand side of the 
above relation can never 
for this last figure are 



_7 
above relation can never be made less than about 10 . Empirical values 





Case I (n 


- 5) 


Eigenvalues s 


1 x 10~ 7 




Eigenvectors? 


2 x 10"' 




Eigenvalues: 
Eigenvectors i 


Case I (n 
3 x 10~ 7 
8 x 10" 7 


= 18) 



Case II (n ■ 5) 
U x 10" 7 
5 x 10" 7 



For an example of the determination of c, suppose that 6 figures were 
designed in the eigenvectors and that |a - Xj 1Bax xv 3« Then 



1 ~7 

i x 10 ' 



2 C X3 x — r- - 150 x 10~ 9 C256 x 10~ 9 =2 8 2~ 30 



2 



-S4 



and therefore c = -24* (A useful relation is that 2 * 1024 " 10**). 
max 

In practice , one would probably use +c ■ -SO to be certain of the desired 
accuracy o As is implied , the value of c can always be safely lowered 
(except not below approximately -60) since the only effect will be to 
raise the computation time somewhat (i e at worst up to 50 or 10Q£ more 
time). For case II when the eigenvalues of S and H are both of the same 
order 9 c^ and c„ should be about the same, since the overall accuracy 
is determined by the least accurate diagonal izationo 



Form of Output 

All of the program output is given photographically with each 
secular equation solution beginning on a new frame „ The various displays 
may be divided into three classes s 1) the initial data, 2) the inter- 
mediate results (case II only) 9 and 3) the final results. Each display 
of the jfirst two classes will consist of a "heading", line followed by 
either a "symmetric" or "square" pattern of the pertinent numbers* The 
third class consists of a heading line followed by the "rectangular" 
pattern of X, in the first column and then the *y"«i *° "Vt n ' eac ^ * n a 
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column, starting from the third column. The first class is normally dis- 
played but may be suppressed , the second class is normally suppressed but 
may be displayed, and the third class is always displayed. The program 
will use as many frames as are necessary to include all the displays, 
starting on a new frame if a particular display will not entirely fit on 
the remainder of a frame. 

The "heading" line will have in the first column, the identifying 
code number of the solution, kf in the second column the criterion value, 
c„ (except on the display of S)§ in the third column the criterion value 
c-, (except on the display of H) § and in the fourth column, an intermediate 
result identification number, i, (only for class 2 displays). Case I 
results may be distinguished from Case II results by the absence of c_ in 
the result heading line c 

The selection or suppression of the various displays depends on 
the contents of the associated suppressor registers; +0 for display! 
-0 for suppression. These displays and their suppressor registers ares 





I 


display 


Contents 


Pattern 


Decimal Address of 










Suppressor Register 


Class 1 


{ H 


H 


sym. ] 


159S 




Is 


S 


sym. J 






f 1 


S f 


sym. 


1599 




u 


u s» 


sq. 


1600 


Class 2 


y 


S-V2 


syaio 


1601 






k 


HS" 1 ' 2 


sq. 


1602 






5 


H' 


sym. 


1603 






,6 


V 


sq. 


1604 


Class 3 


{ 


x+y 


Vf 


rect. 





Symbolically, we have for the various forms of output: 



"Heading" 



"Symmetric" 























• — 






k 




C H 




C S 




i 





























8 



(continued on next page) 



ra-58 



Page 8 of 12 pages 



"Symmetric tt 
(continued) 



"Square 1 * 



"Rectangular" 



11 



12 



13 



U 



15 



9 



10 



1 




n+1 




2n+l 




3n+l 


o o s 

n 




o o • 

2n 




» o o 

3n 




• a • 

4n 
















4n+l 


• 


o o o 

5n 




* • 


e 












1 

o o o 

n 




1 

o e » 

n 




n+1 

e • o 

2n 
















2n+l 

o o o 

3n 




3n+l 

o o e 

4n 


o o o 



The total number of film frames required per solution depends 
on the number of displays not suppressed but may be computed from the 
schedules given belowo If two or more displays are to fit on the same 
frame, they will be separated by two additional spacing lines beyond that 
given in the schedules {For n>S, no two displays will fit on the same 
frame.) Each frame can contain a maximum of 36 lines* 



matrix order 
_ 

2 
3 



number of lines (or frames * f) 



heado + sym 

3 

U 

5 

(continued 



head 



sq, 



3 

K 

5 

on next page) 



heado + rect. 

4 
9 
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A 


6 


6 


11 


5 


12 


13 


13 


6 


13 


15 


15 


7 


U 


17 


25 


8 


15 


19 


28 


9-10 


f 


f 


f 


11-12 


f 


2f 


2f 


13-U 


2f 


2f 


2f 


15-16 


2f 


2f 


3f 


17 


2f 


3f 


3f 


18 


2f 


5f 


5f 


19 


2f 


5f 


6f 


20 


3f 


5f 


6f 


21-22 


4f 


6f 


6f 


23-24 


4f 


6f 


7f 


25-26 


• 5f 


7f 


7f 


27-28 


5f 


7f 


8f 


29-30 


6f 


8f 


8f 


31-32 


6f 


8f 


9f 



For convenience the most common situations are summarized: 





number 


of frames 


per solution 


matrix order 


normal 
case I 


normal 
case II 


Pull suppression of 
cases I and II 


1-4 


1 




1 




1 


5-6 


1 




2 




1 


7-8 


2 




2 




1 


9-10 


2 




3 




1 


11=12 


3 




4 




2 


13-U 


4 




6 




2 


15-17 


5 




7 




3 


18 


7 




9 




5 


19 


8 




10 




6 


20 


9 




12 




6 



(continued on next page) 
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21-22 
23-24 
25-26 
27-28 
29-30 
31-32 



10 
11 
12 
13 
U 
15 



H 
15 
17 
18 
20 
21 
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6 
7 
7 
8 
8 
9 



Information for Filling Out Performance Requests 

1* Time estimates: These are difficult to state in general but 

3 1 
roughly T TT ^ 2T T and T T = fir where K^rr sec. Variations 
^ * case II case I case I 10 

from this formula depend on the difficulty of diagonalization (i.e. on the 

size of the off -diagonal elements and the stringency of the criteria).. * 

2o Program stops Stops automatically on siO in register 41 
octal. 

3 » Camera output; See film frame schedules above for the 
number of frames per solution o Camera first used in see- 

4o Magnetic drum: Case I uses auxiliary drum groups 1 9 2, 3$ 
case II uses auxiliary drum groups 1, 2, 3> 4-« Drum first used in sec 

5» Operating instructions: These are for a single-parameter 
tape 

E* fb A 9 RI, 
fb 10172-20-B, RI, 



f c C 

where A ■ Mummy logging tape tt § 

(tape room will prepare 
this) 



9 RI> RS« 

C 334? (case I) 
1 331, (case II) 

G * parameter tape number 



B 



For multiple-parameter tapes 9 the corresponding instructions are 

E, fb A , RI 9 
fb 1072' jj 20«hB^ RI> 
Turn off si 1 switch 
fc G 9 RI, 

* An additional time of about 10 seeonds must be added for each full 
frame displayed. 
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Program Alarms 












Probable 




Type of Alarm 






Octal Identification 


Immediate Cause 


Cause 


1. 


Check register 








PC « 3764 
AR « 3453 


GD no. unscale- 
factored 


a 


2. 


Check register 








PC ■ 3764 
AR » 3503 


GD no. in MIA too 
large to store 


a,b 


3o 


Check register (case 
only) 


II 




PC - 2601 


An eigenvalue of 
S is<0 


a,b,c 


4-o 


Divide error (case 


II 


only) 


PC - 3612 


An eigenvalue of 
S « 


a,b,c 



Probable Causes 

a. GD number of initial data out of phase, due to either a missing decimal 
point, a missing seventh hole in the tape, an incorrect storage address, 
or a faulty parameter tape read- in by the computer. Compare print of 
data tape with display of input data. 

b. Input data is very poorly scale-factored or, (case II only), one of 
the eigenvalues of S is extremely small. 

Co The matrix S is not positive definite* 



Additional Remarks for Experienced Programmers 

1. If matrix elements are to be generated, decimal registers 
38 to 1563 inclusive are available for programs and will be restored after 
each solution. Additional storage is available on the auxiliary drum 
starting at the decimal address 7630 (case I) and 10206 (case II). All 
generation program tapes must have a NOT PA included; the PA already in 
the program contains buffers b and 6b inclusive and a single cycle counter. 
Thus all cycle instructions may be used, but no "isc" orders (except isc0)< 
Buffers b through 3b are used as temporary storage during a solution and 
4b through 6b are unused. 

2o The program also contains a DIB/DOB subroutine 5 the decimal 
entries arei 
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for DIB 




for DOB 


spl605 




spl6l0 


[cm] 


core memory address 


[cm] 


M 


auxo drum address 


[da] 


w 


no. of rego to transfer 


w 


< 


- return point in wi mode — 


— ¥■ 



lyuaf 



F.J. Corbato 
March 15, 1955 



